clear; %Deletes all variables
clc; %Clears the log screen
T = 10; %Number of observations per simulation
K = 10; %Number of resamplings per simulation
N = 1000; %Number of simulations
% For simplicity, we set var(epsilon) = var(x) = 1
betas = NaN(4,N); %Empty matrix to store betas
for n = 1:N
 % CHOOSE ONE OF THE FOLLOWING DISTRIBUTIONS. COMMENT OUT THE OTHER
 % I'M SIMULATING X AND EPS TOGETHER FOR SIMPLICITY
 dat = randn(T,2); %Normal distribution
%  dat = trnd(4,T,2); %Student's t-distribution with 4 dof
 x = dat(:,1); eps = dat(:,2); %extract x and eps
 % CHOOSE ONE OF THE FOLLOWING FUNCTIONAL FORMS. COMMENT OUT OTHER
 % y = x + eps; %true model for y is linear
 y = x.^2 + eps; %true model for y is quadratic
 
 rhs = [ones(T,1) x]; %set up rhs matrix for ols
 betaols = (rhs'*rhs)\(rhs'*y); %estimate via OLS
 % DO THE RESAMPLING PROCEDURE
 betas_tmp = NaN(2,K); % Empty matrix to store betas
 for k = 1:K
 dat_tmp = datasample([x y],T); %resample with replacement
 x_tmp = dat_tmp(:,1); y_tmp = dat_tmp(:,2);
 rhs = [ones(T,1) x_tmp]; %set up rhs matrix for ols
 beta_tmp = (rhs'*rhs)\(rhs'*y_tmp); %estimate via OLS
 betas_tmp(:,k) = beta_tmp; %write result to matrix
 end
 betas_tmp = mean(betas_tmp,2); %average betas across resamplings
 betas(:,n) = [betaols; betas_tmp]; %write ols and resampled betas to matrix
end
mean_betas = mean(betas);
var_betas = var(betas);
sk = skewness(betas);
ku = kurtosis(betas);